      SUBROUTINE  INTEG
      REAL*8  RAD,TMASS,GS,RHO0,GA
      REAL*8  D,H,RHO,ELAMB,EMU,GR,XI,PHI,ETA
      COMMON/MODEL /NLAY,ISO,LAME,NDIV,DMAX,ILAY,
     1              RAD,TMASS,GS,RHO0,GA,
     2              HI(500),RHOI(500),VPI(500),VSI(500),NAME(20),
     3              XII(500),PHII(500),ETAI(500),
     4              D(1000),H(1000),RHO(1000),ELAMB(1000),EMU(1000),
     5              GR(1000),XI(1000),PHI(1000),ETA(1000)
      REAL*8  RMAX,RC,HC,H1,H2,ALP,WN,WN2,C,C2,FRQ,FRQ2,T,U,ENGY,DELTA
      REAL*8  A,F,Y
      COMMON/VALUE /MODE,IERROR,I,ISTEP,J,K,L,N,NSOL,ISUM,KMAX,
     1              NMAX,IBOTM,ITOP,LY,LD,ACR,ELLIP,
     2              RMAX,RC,HC,H1,H2,ALP,WN,WN2,C,C2,FRQ,FRQ2,T,
     3              U,ENGY,DELTA,
     4              A(6,6),F(20),Y(6)
      REAL*8  YN,YB,SUM,Q
      COMMON/SOL   /YN(6,1000),YB(6,3,20),SUM(20),Q(3,21),
     1              WNB(100),TT(100),CC(100),UU(100),ENG(100),ELL(100),
     2              AC(100)
      REAL*8  Y1,Y2,Y3
      COMMON/Y123  /Y1(6,3),Y2(6,3),Y3(6,3)
      REAL*8  W,WW
C
C PURPOSE - TO INTEGRATE DIFFERENTIAL EQUATION UPWARD
C INPUT - J=-1 WHEN SEARCHING FOR AN EIGENVALUE
C          = 2  WHEN CALCULATING EIGENFUNCTIONS
C         ITOP,IBOTM= UPPER AND LOWER LIMIT OF I
C
C INITIALIZATION
C
      I=IBOTM
C
      CALL  INIVAL
C
      IF(IERROR)  900,1,900
    1 CONTINUE
C
      DO  2  NS=1,NSOL
      DO  2  LL=1,L
      Y1(LL,NS)=YB(LL,NS,1)
    2 CONTINUE
C
      IBOTM=I
      NN=0
      N=NMAX+1
      K=1
      ISTEP=-1
C
C
  100 CONTINUE
      NN=NN+1
      N=N-1
      IF(J-1)  200,200,110
  110 CONTINUE
C
      CALL  EIGFUN(Y1)
C
  200 CONTINUE
      IF(I-ITOP)  1000,1000,210
  210 CONTINUE
      IF(H(I))  500,300,500
C
C DISCONTINUITY
C
  300 CONTINUE
      I=I-1
      IF(EMU(I+1))  310,320,310
  310 CONTINUE
      IF(EMU(I))  100,330,100
  320 CONTINUE
      IF(EMU(I))  350,100,350
  330 CONTINUE
      GO TO  (340,350,340,350,340,350),MODE
C
C LOVE WAVE TYPE
C
  340 CONTINUE
      I=I+1
      ITOP=I
      GO TO  1000
C
C DISCONTINUITY OF DIFFERENTIAL EQUATION
C
  350 CONTINUE
      K=K+1
C
      CALL  CONT
C
      GO TO  100
C
C NORMAL SEQUENCE - TO INTEGRATE 2-STEPS
C
  500 CONTINUE
C
      CALL  STEP
C
      IF(J-1)  700,700,600
C
C ENERGY INTEGRAL
C
  600 CONTINUE
      N=N-1
      CALL  EIGFUN(Y2)
      N=N-1
      CALL  EIGFUN(Y3)
      N=N+2
      CALL  ENERGY
      N=N-2
C
  700 CONTINUE
      I=I-4
      NN=NN+2
C
      DO  710  NS=1,NSOL
      DO  710  LL=1,L
      Y1(LL,NS)=Y3(LL,NS)
  710 CONTINUE
      GO TO  200
C
C
 1000 CONTINUE
      NMAX=NN
      K=K+1
      KMAX=K
      IF(J-1)  1100,1100,1200
C
C STORE SURFACE VALUES
C
 1100 CONTINUE
      DO  1110  NS=1,NSOL
      DO  1110  LL=1,L
      YB(LL,NS,K)=Y1(LL,NS)
 1110 CONTINUE
C
      IF(J)  1120,900,900
C
C COMPUTE CHARACTERISTIC EQUATION
C
 1120 CONTINUE
C
      CALL  CHAREQ
      GO TO  900
C
C NORMALIZATION OF EIGENFUNCTION
C
 1200 CONTINUE
      W=1.0D+0
      IF(MODE.GE.3)  W=RAD-D(ITOP)
      W=W/YN(1,1)
      WW=W*W
      DO  1210  NNN=1,NMAX
      DO  1210  LL=1,6
      YN(LL,NNN)=W*YN(LL,NNN)
 1210 CONTINUE
      DO  1220  IS=1,ISUM
      SUM(IS)=WW*SUM(IS)
 1220 CONTINUE
C
C ENERGY INTEGRALS IN THE HOMOGENEOUS SPHERE
C
      CALL  DENERG
C
      GO TO  900
C
C EXIT
C
  900 CONTINUE
      RETURN
      END
